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We extend a previous statistical mechanical treatment of the traveling salesman problem by 
defining a discrete "site disordered" problem in which fluctuations about saddle points can be 
computed. The results clarify the basis of our original treatment, and illuminate but do not resolve 
the difficulties of taking the zero temperature limit to obtain minimal path lengths. 
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I. INTRODUCTION 



^1" 

Theories describing disordered materials are based either on disorder arising from randomness in the strength 
of bonds or from randomness in site location. Although most theoretical progress has been made in the bond- 
disordered case, site-disorder is probably closer to the picture of real materials. A similar situation exists in the best 
known problem of combinatorial optimization: the traveling salesman problem (TSP), here taken to be the stochastic 
I . version. The Euclidean problem is defined in terms of randomly assigned city locations, but analytic progress is 
2h ' most developed for the random-link model [J 0| m which inter-city distances are chosen independently from some 
probability distribution. Recently we have presented a functional statistical mechanical technique to compute the 
average length of the TSP path at finite temperature for the site-disordered model [3, 0| • In this paper we define a 
discrete model that provides a more rigorous foundation for our functional approach in the continuum and allows us 
to compute fluctuation terms that illuminate the difficulties of taking the zero temperature limit to obtain minimal 
path lengths. 

The site disordered TSP is a straightforward generalization of the ordinary TSP that populates a graph consisting 
of M nodes or sites with N cities. The graph can be abstract or for geometric applications we could take it to be a 
regular lattice. Site disorder arises since the number of cities on each site, which may be 0, 1 or greater than 1, is 
part of the specification of the problem, and we will have in mind a quenched random distribution. This problem is 
defined in the first section below and various regimes depending on the average density of cities, N/M, are identified. 
■ Section [TTTl shows how the results of our previous functional formalism can be obtained from the N ^> M high density 
^sD [ or "coarse" regime of the discrete site disordered TSP for a fixed rather than random distribution of cities. These 
• results correspond to the leading order saddle point approach, and we show how they can be extended by analyzing 
1 fluctuations in the discrete model. 

Our interest in the fluctuations arises from the possibility that they may be used to probe minimal TSP path 
lengths. This was not possible in our previous work since the leading saddle point term always gives path lengths 
that grow proportional to N and correspond to strictly finite temperature. In certain cases the fluctuation terms have 
the correct scaling and may be expected to dominate in the low temperature regime where the minimal path lengths 
can be accessed. Unfortunately this program is not realized. The fluctuations can be computed explicitly for some 
graphs with simple structure and we discover that in the limits needed to find minimal path lengths the saddle point 
approximation breaks down in a complicated manner. The manner of this breakdown is most clearly exposed for the 
case M — 2 that is analyzed in detail in the appendix. 

The remaining sections of the paper discuss the site disordered problem in two other regimes. Firstly the case 
N = M, precisely the ordinary TSP, where it is interesting to see that the saddle point approximation breaks down 
in a completely different way from in the N 3> M coarse regime. Secondly the case N <C M, which for a regular, 
dilutely occupied lattice can also correspond to the ordinary geometric TSP, and we discuss the significant interest of 
this regime, and point out the additional techniques that may be needed to resolve the problem in this regime. 



II. FORMALISM 



Consider a graph consisting of M sites with a matrix Jjj specifying the distance between any pair of sites i,j = 
0, 1,2 ... M — 1. For simplicity the diagonal terms 3 a — Jd are taken all equal but the matrix need neither be 
symmetric nor positive though we shall generally take it to be so. Imagine that there are rij cities {rij > 0) at site 
j. We will call these quantities the "occupation numbers" and for a total of A^ cities they obey X^jlo 1 n j = N. A 
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path through the N cities is specified by the sequence of site indices jo,ji ■ ■ -Jn-i, so depending on the occupation 
numbers (n/s), some sites might never be visited, or visited more than once. A particular path, jo,ji,j2 ■ ■ -Jn-i, 
can also be characterized by an element of the permutation group, <r £ Ejv, consisting of a permutations of M site 
indices, each repeated rij times, so as to have a total of N indices. In this notation the path is written <7o, o\ . . . <Jn— ii 
and the cost or total length for this path is: 

JV-l N-l 

D(a) = Y, J-r m ,a m+1 = E J ^ m+1 ■ (1) 
m— m— 

Where we identify jjy with j . As is often the case in this kind of problem, the cost is independent of cyclic changes 
and also (if is symmetric) reversals of the permutation. 

A full statement of the problem involves a specification of the occupation numbers as well as the cost matrix between 
all sites. Because we take a site-disordered approach we shall generally assume that the cost matrix is a fixed known 
quantity and will often have a geometric basis in mind for which the cost matrix follows from a regular lattice graph. 
Disorder can appear through the occupation numbers on each site. 

By varying N/M, the average number of cities per site, we see that this model describes various regimes. If we 
take N = M and set all the = 1, then we immediately return to the TSP as usually formulated with arbitrary 
cost matrix that may, or may not, have a geometric basis. When we consider geometric problems such as the 
Euclidean TSP it is numerically advantageous to imagine the sites as being on a regular lattice, thus allowing the 
use of integer arithmetic to compute [S] , and by considering sufficiently dilute occupation N -C M the continuum 
result is recovered. In this Euclidean lattice case, the N = M situation is trivial as it corresponds to a fully occupied 
lattice where the optimum TSP paths only involve nearest neighbour links. The approach to the continuum and 
standard TSP as lattice occupation is made progressively more dilute, N <C M, has been considered numerically 
in ((J. Analytically, the most tractable regime is when N ^> M and this may be regarded as a coarse, in the 
renormalisation group sense, approximation to the TSP in which groups of close cites are bunched together at the 
same site ignoring the detailed distances between them. Even in this coarse regime, information about the standard 
TSP can be obtained by arranging for the diagonal distances Ja to be small. 

The statistical mechanical partition function for this problem at temperature 1//3 is defined by the following sum 
over permutations that visit sites the correct number of times. 

ZnW, {%•})= £ exp (-/3.0(a)) (2) 

This partition function Zn((3, {n,j}) depends both on occupation numbers and Jy's and therefore encompasses both 
bond and site disorder. The sum over permutations may be written as an unconstrained sum over all possible paths 
of N links with a Kronecker delta constraint to ensure that site j is visited precisely rij times. 

M-l M-l / N-l \ 

Z N (P,{nj}) = £ II n * 1 S n k ,Zl-J s k , ]m ex P -P E J ^ m+1 • (3) 

iotil---JJV-l=0 k=0 \ m=0 I 

The constraints appearing here should not be confused with the constraints that forbid disconnected paths in the linear 
programming approach to the TSP; here the path is always connected, but must be forced to visit sites the correct 
number of times. The combinatorial factors of n^! arise from the different possible orderings of visiting multiple cities 
at the same site, and are necessary to ensure that with no weights, (3 = 0, we recover Zn — N\, corresponding to the 
total number of possible orderings of N cities. It is implicit in the formalism that "Y^n^ = N since this constraint is 
enforced by one of the Kronecker deltas. 

Writing each Kronecker delta in terms of an integral representation we find, 

M-l r M-l j_ t _n*-l . JV-l /N-l \ 

w k» = e f n k 2m n jr ex p h 3 £ j — ( 4 ) 

3o,j'l--iw-l=0 k=0 " m=0 jm \ m=0 J 

M-l „k M-l ^ ^ , / M-l N-l \ / JV-l \ 

E / II ~kr 1 exp r E nk '^ k ' ~ i E ) exp ^ E J M m +i ) ■ ( 5 ) 

io,h-jN-\=0 * k=Q \ fc'=0 m=0 / \ m=0 / 



Where the second form is obtained from the contour integral on the first line by setting Zk — e l ^ k . Although the 
integrals are finally evaluated using this second form, it is helpful to keep the first form in mind to check what 
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deformations of the contour are possible in the presence of branch cuts or poles. Rearranging, the partition function 
can be written in terms of an integral over a subsidiary problem, 

Z N (f3,{nj}) = I J[ k k ' exp i ^ n k ,^ k , Z N (/3,{nj}) (6) 

J- 71 k=0 11 \ fe'=0 / 

The subsidiary statistical mechanical problem does not depend on the occupation numbers, n k , and is defined as 

M-i / N-l N-i \ 

z N (p, {ti,}) = ex p (- 1 E - $ E J ^^ +1 = Tr tN - w 

jo ,ji-..jjv-i=0 \ m=0 m=0 / 

Where the M x M, transfer matrix, T is the symmetrized form 

T = RBR. (8) 

R is the diagonal matrix with diagonal elements R k — z k = exp(— i/^/2), and B = cxp(— /3J) is defined in terms 
of the distance matrix. In this paper we shall only consider symmetric distance functions, so the transfer matrix is 
symmetric, though generally complex and therefore not hermitian. 

Distance matrices that have some kind of regular geometric basis yield additional structure for B. In one dimension 
this is most obvious, for example a regular one dimensional lattice with J^- cx makes B Toeplitz, and if the lattice 

has periodic boundary conditions B is also circulant. For regular lattices in two dimensions B is block circulant with 
each of the blocks Toeplitz, or themselves circulant in the case of periodic boundary conditions. In higher dimensions 
there is a hierarchical structure of B block circulant with blocks that are themselves block circulant. For all these 
regular lattices the eigenvectors are related to discrete Fourier transforms 0] We also consider fully connected M- 
dimcnsional graphs with = Jjy + (Jd ~ JN)$ij, which are particularly tractable. We often allow a diagonal term 
in the cost matrix, as this can be used to probe certain aspects of the theory. We will not consider asymmetric cases 
in this paper, but examples such as the wallpaper problem [8] with = 9(i — j) can be treated. 

The solution of the subsidiary problem can be written in terms of the eigenvalues A a of T, 



M-l 

Z N (P,{ H })=TrT N = J2 X p ( 9 ) 

p=0 



For general values of the integration variables z k or /i k , the eigenvalues will be complex. Under the scaling z k — > az k , 
eigenvalues transform as A — > A/a, and below we shall find that this scaling symmetry is related to the Lagrange 
multiplier mode enforcing the constraint ^ n k — N . 

To close this section on the formalism, let us interpret it in the case of the N — M regime with all occupation 
numbers n k = 1, corresponding to the usual TSP. The partition function then becomes, 

Z N (0)=<fl[pTr(T* f ) (10) 

With T = RBR and R 2 = diag(zQ , z^ 1 . . . z^ 1 ^), the trace is a polynomial in the variables zj 1 that acts as a 
generating function for paths that visit site j a number of times given by the power of zj . The weight associated 
with a particular path is the exponential of the length of the path. Since all paths visit ./V sites, every term in 
the generating function polynomial has equal degree N . Depending on the symmetry of the distance matrix, the 
generating function often has some symmetry in the variables Zj. The integration defining the partition function for 
the TSP picks out the pole corresponding to the term Y[ zj 1 m the generating function, thereby requiring each site 
to be visited exactly once. From this point of view, the approach is similar to other combinatorial problems where 
the generating function is known, and in many graph theoretic problems a similar integration can be used to select 
the desired term. For example in the mathematical literature, the saddle point technique has been used to obtain 
asymptotic formulae for the number of Eulerian oriented graphs 0. For (3 — 0, distances are irrelevant and vertices 
are either connected or not, so the system becomes purely topological and the problem is to count Hamiltonian paths 
for a graph with incidence matrix B. The formalism also resembles the field theory for studying random Euclidean 
matrices 1 101 - 
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III. COARSE TSP REGIME: N > M 



Consider the limit of many cities, N — > oo, with M fixed. In this limit the occupancies, n^, scale with N, and it is 
convenient to rewrite them in terms of a "density of cities" , pk where: 

N 

n k = j^Pk- (11) 

The normalization of pk is ^ fe=0 Pk = Ad ■ 

In a geometric context, the formalism is now based on densities of cities rather than their explicit locations and 
this limit might be regarded as viewing the system at coarse resolution, smearing out the precise city locations. Such 
coarsening techniques are indeed used in heuristics for solving TSP instances, although in a different way, for example 
the Held Korp heuristic patches together paths from different cells. 

In this regime a diagonal term Ja = Jd in the cost matrix can have a special significance. If this self connection 
corresponds to a distance larger than typical neighbour distances, then it will have little effect at low temperature 
and the optimum path for uniformly distributed cities may correspond to the optimum of the standard TSP, repeated 
N/M times. If on the other hand, Jd is short in comparison with other distances, then the favoured paths at low 
temperature will satisfy the requirement to visit cities n k times by simply looping that many times before proceeding 
to a neighbour. In this case the true TSP path is found simply by ignoring the contribution from the loops. For 
various simple models with small M, we have considered this small Jd situation explicitly. At low temperature, by 
regarding the neighbour couplings as corrections to the diagonal part of B then a perturbative approach can be used 
to compute the eigenvalues of T and hence the contribution to the partition function which correctly predicts the true 
TSP path. The hope of being able to compute such terms for arbitrary M is a motive for studying the fluctuations 
about the saddle point. Unfortunately, it appears that the saddle point method, which is our main tool, generally 
fails to give the correct subleading term in this limit for technical reasons that are discussed in the appendix and in 
the sections where particular models are considered. 

By inserting the solution of the subsidiary problem into the expression for the full partition function, we find, 

/M-l \ M-l 

Z N (f3,{ nj }) = exp ( J2 lo S^ ! J E Z N^,{pj}) (12) 

M-l , I . M-l \ 

fc=0 \ fc'=0 / 

Where the contributions from each eigenvalue are now distinguished, at least by a label p = 0, 1, 2 . . . M — 1. In 
practice, globally labeling the eigenvalues of T as the zj. vary is not straightforward. 

By writing in terms of pk the N scaling in the exponent is exposed that allows the integrals for each Z p to be 
evaluated using a saddle point approximation For the simplest application of this technique to be valid, the 

eigenvalue must have no N dependence so we are forbidden to scale the temperature with N [ijj]. The saddle point 
equations for each Z p are, 




1 d\ p ; p k 



( 14 ) 



A p dp k M 

Now, by taking derivatives of the eigenequation, Tv p = A p v p , noting that dT^/Opk = —i{$ki + <5fcj / 2 , recognizing 
that the left and right eigenvectors are the same due to the symmetry of the distance matrix and using normalization, 
v p .v p = J2k v k v k = 1' we fi nc l l ne general relation, 

^ = -« 2 (15) 

So the saddle point equation may be rewritten as, 

K) 2 = %■ (16) 

The right hand side of this saddle point equation does not depend on the eigenvalue index p. Ostensibly we may 
choose signs independently for each component of the eigenvector to generate distinct solutions for each Z p , but we 
shall discover that the role of the signs is not transparent. The eigenvector is written vf. = £f Wpk/M, where the 
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signs = ±1. Now decomposing the transfer matrix and writing = £P e l /W 2 \J pkX/M 7 the signs disappear and 
the eigenequations become, 

S ! = iVe^4. (17) 



For all reasonable choices of the distance matrix this equation has a solution with real positive Sk (to denote this 
solution we omit the index p), that can be found by a numerically stable iteration procedure. In many of the explicit 
models we consider, this solution is in fact independent of the index k; for example, in the case of circulant J and 
Pk = 1, this positive solution is s k = s;s 2 = (1/M) i-B-y- For certain particular choices of the distance matrix, 
additional solutions to (jTTJ) may be found. For all cases considered, these additional solutions, while not spurious, are 
less important in the large N limit than the real positive solution noted above. These solutions are discussed in the 
sections below that deal with particular models. 

In the remainder of this section we shall only be concerned with the real positive solution s&. The location of the 
saddle points for this solution in terms of integration variables is given by: 

s 2 

Xz k = M— (18) 
Pk 

Note that all reference to the signs (1 has disappeared. Evidently, since the right hand side is real and positive, there 
is a saddle point on the real Zk axis with Re pk = 0. At this saddle point the transfer matrix is real and symmetric 
and has real eigenvalues. Moreover (provided that the cost matrix is positive) the Perron Frobenius theorem assures 
us that the eigenvector v k = +yj Pk/M, where £j£ = +1, corresponds to the maximal eigenvalue. 

Notice that under the scaling symmetry Zf. — ► az k , eigenvalues transform as A — ► A/a, so the combination Xz k 
appearing in p8[) is invariant. As mentioned above, this scaling corresponds to the Lagrange multiplier mode involving 
a simultaneous shift of all pk, that enforces the constraint ^rik = N. Since this constraint is implicitly satisfied by 
the pk that we specify, there is no need to integrate over it, and we may choose to set ^ \i k = or ]1 z k = 1. As 
we shall see explicitly in section IV, the stationarity condition associated with this mode is not necessary to solve the 
saddle point equations since the normalization condition for p k can replace this information. For the real solution, s/-, 
the requirement Y[ z k = 1 forces A M to be real. We can therefore identify a total of M saddle points simply related 
to the real one by overall phases. These saddle points, arising from the positive real solution to (fT7| will be termed 
dominant. 

z (p) = e 2^/M z (o) ( p = 0, 1 . . . M — 1 (19) 

.(0) 



4 ■ mr 



tlt-\ 2 \ 1 /M 
Si 



A<°> = M TT — (22) 
\t=o Pk) 

Note that for the case in which p k — 1 and the solution s k is independent of fc, these equations simplify dramatically: 
\(°) — Ms 2 and the saddle points are equally spaced around the unit circle. 

It is tempting to identify each of these M saddle point solutions with a single dominant saddle point for each of the 
M integrals Z p . In that case, the phase factors cancel in the leading contribution and each Z p contributes equally to 
the total partition function, resulting in an overall factor of M that can be interpreted as arising from the redundancy 
of cyclic permutations or equivalently the choice of origin of the paths. However, the derivation above does not 
provide an unambiguous identification of the saddle points with the Z p since the eigenvalues are not globally labeled. 
A suitable approach would be to label eigenvalues according to their order at some reference point, say z k = 1 where 
they are real, and then to smoothly deform zj~. This procedure is subject to difficulties since the constraint Y[ z k = 1 
should be maintained whilst deforming, and incorrect conclusions arise if one takes all Zk = e lB and varies 9 = to 
9 = 2-np/M. Some intuition can be gained from explicitly considering small M examples and in the appendix the 
TSP on two points is treated. In that case, although the picture is complicated by the existence of another solution 
to (|17[) . the identification is correct, with the dominant contribution (from the positive real solution) yielding a single 
saddle point for each of the two eigenvalues. 

Henceforth we shall assume that the picture of saddle points given above is correct. In the next subsection we 
consider the order 1/N fluctuations about the saddle point, but first we deal with the leading contribution from the 
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dominant saddle points. The free energy is, 



N 



(3F(f3, {p k }) = log Z N = — Pk log \z k + TV log 



TV log TV 



(23) 



Where we have used the Stirling approximation for n^!. 

The expected path length, (E), and other quantities such as correlations along the path can be computed from 
the partition function by appropriate derivatives. The expression for the path length can be rewritten in terms of s k 
itself, but higher moments explicitly involve derivatives. The final equations determining the leading contribution in 
this coarse regime are: 



if 2^ e 



AI 



[3F = -\ogZ N = -27V-L^p fc log Sfe -TVlogTV 



(E) 



TV 
M2 



PkPk'Jkk' e ^ 0Jkk , 



kk> 



SkSk' 



(24) 
(25) 
(26) 



Notice that the expected path length is always order TV in the limit we have considered. This is indeed expected at 
finite temperatures where the number of relevant contributing paths grows as a factorial of TV. To see this we may 
compute the entropy as: 



— ^Pfclogs fc 

k 



NlogN = NS E + NlogN 



(27) 



Where in the final form we have isolated the extensive, Se, from the non-extensive contributions to the entropy. Then 
by inverting the expression for the path length obtained from equations (|26[) . the entropy can be re-expressed as a 
function of path length. The significance of the entropy in this problem is that the number of paths of given length 



Af(E) 



TV! e 



NS E 



(28) 



The factorial growth in the number of paths with the number of cities characterizes this TV 3> M regime and indicates 
that this is not the regime where the optimum path can be found. To approach that regime, where cities are linked to 
one of their nearest neighbours thus making the number of paths grow exponentially, some scaling is necessary. Either 
the distances Jy can be scaled so as to preserve constant density of cities as TV — > oo, or equivalently the temperature 
can be scaled with TV pL3[ . This kind of scaling would invalidate the simple use of the saddle point to solve for the 
partition function as the large parameter TV would not appear as a simple factor in the exponent. Therefore the 
leading order of this coarse model cannot access the optimal paths of interest in the TSP. However, as we recall below, 
the approach can be used at finite temperature and in studying the maximal TSP problem [14j | where no scaling is 
needed to find the regime of maximal optima. 



A. Continuum Geometric Models 



In order to relate the leading order site disordered TSP results above to the results of our previous functional 
approach we take a continuum limit with indices j corresponding to points Tj in some d dimensional domain. The 
distance matrix is defined by the Euclidean distances between the points, = |r,; — rj\. Assuming a domain with unit 
volume, the appropriate lattice spacing scale is a = M^ 1 ^ and the translations necessary to obtain the continuum 
limit are 




nk ~* T7P^ rk ) ( 30 ) 



Where p(r) is the continuum density of cities. 
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FIG. 1: Theoretical prediction for the average path length E for two dimensional continuum TSP on a box, sphere and torus 
as a function of /3. Negative ft corresponds to the MAX TSP. 



In the simplest procedure, the points are laid down on a regular lattice and the continuum limit is straightforward. 
In this case there is no randomness whatsoever. We could also imagine M points laid down at random on the domain, 
then some weight factor is needed as we take the continuum version of the discrete sum. This weight factor appears 
in exactly the same places that the density of cities p does and in the end they are the same. Irrespective of the 
way the continuum limit is taken, the formalism only depends on the density of cities and all information about the 
precise location of the sites disappears. 

The equations now become, 



s (r) = J d 
f3F = -2N J d d rp(r)logs(r) - NlogN 



d -ff|r-r'| P( r ') 

S (r') 



(E) 



N / d d rd 



d ,p(r)p(r') 



-/3|r-r'| 



(31) 

(32) 



(33) 



s(r)s(r') 

These are the same equations written down in 0] using a functional formalism in the continuum. We now briefly 
summarize some of the results for low dimensional models based on these equations. The most interesting predictions 
of this theory are for the MAX TSP 14] corresponding to negative ft. 

The formalism presented above yields the correct high temperature expansion as a series in ft, and has also been 
tested by comparison with Monte Carlo simulation of the TSP at finite temperature |J|. Except for the case of 
closed symmetric domains to be discussed below, we have not found any non-trivial analytic solutions of the nonlinear 
equations (|33[) . For the special closed domains, such as a disc and torus in two dimensions, a constant solution exists. 
The significance of this observation is that the annealed approximation, in which the average over cities is taken at 
the level of the partition function Eq. ([2]), is exact for these domains. 

For the sphere, the annealed/quenched equations have a constant solution with: 



-in 



d 2 r e-? e = 



2vr 



S-2 



(ft 2 + Air) 



leading to an expression for the average path length, 

E/N 



2ft 



(34) 



(35) 



/3 2 +4vr 2(ev^/3/2 + l) 

Note that in the limit ft — > — oo, the path length per link becomes the half circumference corresponding, as is the case 
for all the closed domains, to the maximum distance two points can be apart. 

For a torus we find that the equations yield 

p e -0/2 < 



frr/4 



2cos6> 



-)d6 



(36) 
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As [3 — > — oo, the integral can be approximated and again the path length per link becomes the maximum distance 
two points can be apart, E/N — > 1/V2- 

For most domains, including the traditional square box domain, the quenched result is different from the annealed 
approximation and there is little hope of a general analytic solution. In this case, our primary tool is the iterative 
numerical solution of Eq. (|33[) which is stable and can be solved to any required accuracy. 

Figure (QJ shows the average path length for each of the three two-dimensional domains we have considered. For 
the sphere this is given by (f35|) . for the torus it is based on numerical integration of (|36|) and for the box we resort 
to an iterative solution of the original quenched equations. The accuracy of this iterative technique is confirmed by 
reproducing the results for the other domains. In all these cases we have also performed Monte Carlo simulations and 
obtain excellent agreement with the theory. At large positive the topology starts to become unimportant and each 
domain has average path length ~ 2 N / (3. In dimension d we expect E ~ dN / f3, so this leading term vanishes at low 
temperature. 



B. Fluctuations 



Despite the inability of the leading term of this approach to access minimal length paths, it is interesting to consider 
the fluctuations which can contain terms with the correct scaling form. Gaussian fluctuations about the saddle point 
yield subleading corrections in A, containing a factor det -1//2 H, where H is the matrix, 

d 2 x dv\ 

ilfefe' = --r ttt = q \SP- (37) 

In this coarse regime, the order of the corrections from these fluctuations will be 1/A with respect to the leading 0(A) 
term in the free energy and will thus correspond to an A-independent term. In one dimension, such an A-independent 
term corresponds to the scaling of the true TSP optimum. To proceed explicitly requires additional knowledge of 
the distance matrix. We shall consider the special case of pk = 1 and a circulant B where a constant solution s 
exists and the eigenvectors are well known to be related to Fourier transforms 0]. This restricted choice of circulant 
matrix still covers various interesting geometrically based possibilities and the slightly different case of block circulant 
matrix is a simple generalization. Rather than use the general determinant formula above, we proceed by explicitly 
expanding about the saddle points. We assume the picture of M dominant saddles and find that besides the phase 
cancellation at leading order, all saddle points contribute equally at the level of fluctuations, and we only need to 
consider fluctuations about the real saddle point, which in this case is at Zk — 1. To expand about this saddle point 
write, Zk = e Vfe , where /ik is small but may in general be complex. Then the transfer matrix is given by: 

T= (l-^R 1 -iR 2 + ...^B (^I-^R 1 -iR 2 + ...^ =T( )+TW+T( 2 )... (38) 

Where the matrix Ri is diagonal and given by: 

Ri = diag(. ..fi k ...) (39) 

and 

T<°) = B (40) 
T (1) = -^(RiB + BRi) (41) 

T (2) = — - (R 2 B + 2RiBRi + BR 2 ) (42) 
8 

The eigenvalue is written in a similar manner as A' -* + A' 1 ' + A' 2 ) + . . . (not to be confused with the labeling of the 
different eigenvalues of T), and we shall usually drop the superscript on the unperturbed A^ ' = Ao = Ms 2 . The first 
and second order corrections can now be found using the standard techniques of quantum mechanical perturbation 
theory (though note that the matrix Ri is not generally Hcrmitian) . It is convenient to use Dirac bracket notation 
for the eigenvectors of the unperturbed T^ ) = B. 

B|fc) = A fc |fc) (43) 

According to the analysis above, the maximal eigenvector |0) has constant components Vj = 1/yM for constant 
density. The choice of a circulant B also fixes the other eigenvectors, to have components vj — e~ 2m: > k / *JM . 
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The first order correction to the eigenvalue vanishes as expected 

m-:i 

M 



AW = (0|T( 1 )|0)=-^^ Mfc = (44) 



fc=0 

At second order, there are two contributions, 



The expression for the partition function (fT5)l involves the following integral over the Gaussian fluctuations 

A I „M-1 

Ztt \ Aq 



Z N = M 



k=0 



Where the terms in the exponent arise from expanding the logarithm and the prefactor M just counts the contribution 
from all the equivalent saddle points. To perform the integrals it is convenient to change from the /ifc's to Fourier 
transform variables Oj in which the exponent will be diagonal. 



M-l 

i 

fJ-k 



M-l 

1= aje 2 ™ jk / M (47) 



3=0 
M-l 

jL V n k e-™M M (48) 
to 

The constant part of /j,k corresponding to the Fourier mode ao imposes the constraint which we have implicitly assumed 
so the integration is omitted, and indeed no ao terms arise in the integrand. Using these variables we find, 



M-l 

(0|T( 2 '|0) = — £ (A + X^au-j (49) 



L x^rx k = m ^ (Ao - Xj) ajaM ~ j (50) 



M-l 



AU - w (xJ-rx-f jaM - j (51) 



It is now straightforward to evaluate the integral over the Gaussian fluctuations since we can take a contour in the 
direction where [i^ is real and therefore a* = aja-j- This corresponds to the imaginary direction in the Zk plane and 
for other dominant saddles the integration should be along the tangent to the unit circle. In this representation the 
fluctuations are already diagonal and it is clear that all the eigenvalues of the fluctuation matrix are all negative. 

/3F = -2Nlogs + -J" log ( Xo + Xj - ) - TVlogiV (52) 

The final fluctuation part is a simple expression in terms of the eigenvalues of the coupling matrix B. Below we 
consider some particular examples with circulant B in this N ^> M regime, the fully connected and geometrically 
motivated problems where expressions for the eigenvalues are available and the computation of the fluctuations can 
be taken further. 
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C. Fully Connected Graph 

Consider the regular fully connected graph Jij = Jjv(l — §ij) + Jo^ij where all sites are connected by equal length 
links and we have allowed a diagonal term. We write e~^ Jo = d and e _/3JjY = n to obtain the circulant form: 



B 



/ d n n 
n d n 
n n d 



n \ 

n 

n 



(53) 



\ n n n ■ ■ ■ d / 

For a uniform distribution of cities pk = 1, the equations (|17p have a uniform solution Sk = s where, 

s 2 = ±(d+(M-l)n). 



(54) 



Saddle points are equally spaced about the unit circle and the eigenvalue \^ = Ms 2 . 

Based on the dominant saddle points associated with this solution, the leading contribution to the expected path 
length is, 



N 



dJ D + (M - l)nJ, 



JV 



d+(M - l)n 



(55) 



At low temperature, depending on the relative size of the diagonal and neighbour distances, we obtain (E) = NJp 
or (E) — NJpj. These results correspond to two regimes in which paths either prefer to revisit the same site or move 
to a neighbouring site. 

As explained at the beginning of this section, we hope that the fluctuation term may yield interesting information 
in the limit Jn > Jd when f3 — > oo, ie n/d — > 0. Fluctuations around the dominant saddle points can be evaluated 
as described above since this problem has a circulant cost matrix. The eigenvalues are very degenerate and besides 
the maximal eigenvalue given above all others take the value d — n. Using the expression (|52[) we find: 



[3F = —N log (d + [M - l)n) + 



M — 1 



log (2d/n + (M — 2)) + (JV 



M - 1 



)logM- JVlogJV 



So in the limit n/d — > the path length is: 

{E) = 



N 



M - 1 



Jd 



M-l 



-Jn 



(56) 



(57) 



In this limit the fluctuations about the saddle point fail to give the expected optimum: {N—M)Jd+MJ^. The reasons 
for the failure concern branch points approaching the saddle points and possibly other (non-dominant) solutions 
becoming important. Indeed other solutions for Sk besides the symmetric solution exist. The details of the failure are 
most clearly exposed in the simplest M = 2 example given in the appendix. 



D. 1-D lattice 



Several models corresponding to a one dimensional lattice with nearest neighbour couplings can be imagined de- 
pending on the boundary conditions and symmetry of the couplings. On a ring with uniform distribution of cities 



Pk = 1, we may simply allow nearest neighbour connections, Ji 
possibility of longer links, Jy = Jjv|i — j\ + JdS^ or = Jn(N - 



JnSz,j+i + Jn^ij-i + JoSij, or consider the 
— j\) + JD^ij, according to which is the closest 
way around the ring. In each case we have allowed a diagonal term, J a = Jjj. The generalizations to problems on a 
line with fixed boundaries do not have a circulant form, but the difference is only a surface effect which corresponds 
to ignoring terms of order e~ M ^/ 2 . 

In the nearest neighbour case we have the circulant form: 



B = 



(d 


n 





n 


d 


n 





n 


d 



\n 



n \ 





d) 



(58) 
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The equations (|17|) then have a uniform solution Sk = s where, 



The leading term in the path length is then: 



s 2 = ±(d + 2n) 



(E) _ dJ D + 2nJ N 
~N~ ~~ d + 2n 



(59) 



(60) 



At low temperature, the picture is similar to the fully connected case. Depending on the relative size of the diagonal 
and neighbour distances, we obtain (E) = NJp or (E) = NJn, corresponding to regimes in which paths either prefer 
to revisit the same site or move to a neighbouring site. 
The eigenvalues at the typical saddle point Zk = 1 are: 



A fe = d + 2n cos(27rfc /M) 

So using (|52p , the free energy may be written: 

Af-l 

2 



Af-l 

(3F = —N log(d + 2n) + -^2 lo S 



fc=i 



d + 2n cos 2 



nk/M\ 
', sin 2 7T k /M J 



+ (N- 1/2) log M- log Nl 



(61) 



(62) 



For a large number of sites, M, the sum can be approximated as: 

Af-l 



fc=i 



d + 2ncos 2 7rA:/M 
nsin 2 7rfc/M 



\ ,w \ A/ /* 2lr /d/n + l + cosw\ 
- Af-l log2+— / dwIogM- 

/ 27T Jo \ 1 — COS UJ J 




= M log 2 + M log - + 1 

n 



The corrections to the path length can then be read off, and in the limit n/d — > we have: 



(63) 
(64) 



E 



N 



M 



(65) 



The finite correction, were it not for the factor of 2, would be exactly what we expect for the ordinary TSP. The 
approach fails again for technical problems with the saddle point approach in this limit as is most clearly exposed in 
the appendix. An additional non-dominant solution for is found in this case also. 

These results can be generalized to one dimensional models with long range interactions and regular periodic lattices 
in higher dimensions with block circulant B. However, in view of the difficulties at the simplest level, we do not pursue 
this route. 



IV. FULLY OCCUPIED TV = M REGIME 

With the occupation numbers set to n& = 1 and N — M we recover the usual formulation of the TSP with a given 
cost matrix Jjj that was discussed in section II. 

.N-l , „k N-l j / N \ 

zn(0) = f n ex p ( n io § a ) = / nf ex p e + n io s ■ ( 6e ) 

J k=0 711 k=0 71 \k' = l ) 

For a regular lattice there is no randomness in this regime and the TSP optimum is a trivial problem. However, it is 
difficult to recover the trivial result from the approach based on the expression above. For example, the 1 dimensional 
lattice can only be solved explicitly for small N. We investigate this case in some detail to illustrate that the saddle 
point approach in principle yields the scaling that would allow a low temperature limit to access minimal path lengths, 
but that other reasons invalidate the approach. Moreover, we hope that results on regular lattices might form the 
basis of an expansion towards more dilute lattices. 
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The saddle point method cannot be applied directly since the factor TV does not multiply all terms in the exponent, 
but consider the following (unitary) change of variables 



N-l 



flic 



a 



3 



\= Y a,e 2 ^ k ' N (67) 

N-l 

J= £ ^ W (68) 



k=0 



This is a 1-dimensional discrete Fourier transform, but note that it is independent of any possible dimension underlying 
the Jy's and holds for arbitrary cost matrices. Unitary normalization has been chosen so the Jacobean is one and we 
obtain, 

Zn(P)= / CX P (zViVao + N log A) . (69) 

Now, by considering the structure of the general transfer matrix, we find that the ao dependence factors out so the 
eigenvalues are of the form A(ao, ai . . . , ajv-i) = e~ ao ^^^ X'(ai . . . , ajv-i)- Therefore all ao dependence cancels in 
the action as we expect since this mode simply imposes the implicit constraint J2 n k — N that we have satisfied by 
setting nfe = 1. Integration over this mode should be omitted to obtain: 

,JV-1 

Zn(P) = J n 2^exp(JVIogA / (oi.02-"OJV-i))- (70) 

A saddle point technique can again be attempted, but we must take care that the eigenvalue itself does not have any 
residual TV dependence. In the case of the TV 3> M regime, this was straightforward as the matrix was M x M, here 
it is ./V x TV, but for local Jij where the number of nearest neighbours does not grow with TV, the eigenvalue can be 
TV independent. 

The stationary conditions from varying at, fc = 1, 2 . . . TV — 1 are, 

N-l 

E e™W N v? = (71) 

3=0 

So, when complemented with the normalization condition on Vk we have, 

vt = 1/TV (72) 

resembling the saddle point condition obtained in the TV ^> M regime. As in the TV ^> M regime there are spurious 
signed solutions of the saddle point equations that do not give a distinct saddle point from the solution with Vk = 
+1/VTV. Proceeding through similar steps to that regime and writing Si = a/A 7 exp(in/2 — \/TVao/2), we obtain the 
final equations, 



N-l 



s k = Y,— (^) 

fe'=0 

0F = -logZ N = -2]Tl°gSfc (74) 



. , T..,p-P J ij> 

(^) = E^— (75) 

33' 33 

Note that no 1/TV factors appear in front of the sums, however, for local cost matrices at sufficiently low temperature 
only neighbours will contribute to the sums, and we find an TV independent expression for Sj. Such a solution 
suggests that the saddle point is valid as it has no untoward TV dependence. Moreover, this constitutes a correct low 
temperature scaling solution since the nearest neighbour distance is fixed as TV increases and the size of the domain 
increases keeping the density of cities constant. We therefore expect that this approach might access sufficiently low 
temperatures to provide information about the optima. 
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To investigate further we take into account the fluctuations about the saddle point, and following the discussion 
for the coarse regime we consider B to be a symmetric circulant matrix. For these special distance matrices the 
saddle point equations have a constant solution z k = z, the transfer matrix becomes T = z B and all eigenvalues 
are directly related to those of B. Using the scaling symmetry of the ao mode we can choose 2 = 1. In fact there 
are a set of iV saddles with z any iV'th root of unity, but they all contribute equally to the integral and we only need 
study the one at z — 1. Determining the Hessian matrix becomes an exercise in second order perturbation theory 
that when aided by the structure of the matrix yields the following quadratic correction to the exponent: 

k=l u K 

Where the eigenvalues A& are of the unperturbed transfer matrix, which at this constant saddle is simply B (Ao is 
the maximal eigenvalue corresponding to the eigenvector with all l's). Note that the zero mode has dropped out 
as expected, and that in the aj~ variables the Hessian is already diagonal. With some care about the direction of 
integration over the saddle, the quadratic correction to the free energy is: 

1 N ~ l /A + A \ 

(3F = - log Z N = -N log A + - E lo § [ ) ( 7? ) 

Note that there are N terms in the eigenvalue sum and typically each term in the sum contributes. For example, for 
the nearest neighbour one dimensional ring the eigenvalues are given by (fBTj) and the sum is order N. The fluctuation 
correction is of the same order as the leading term indicating a failure of the saddle point approach since we can 
anticipate similar magnitude corrections at cubic, quartic etc levels. 



V. DILUTE LATTICE REGIME: N < M 



The dilute regime is the most interesting from the site disorder perspective. We consider a regular lattice with sites 
either occupied or not, n k = 0,1. As the system becomes very dilute, the discrete nature of the underlying lattice 
can be ignored and we recover the traditional TSP. Indeed, numerical TSP optimization codes use integer arithmetic 
following this principle. Here we would be content to investigate the leading correction to the path length for small 
dilution as the numerical simulations Q seem to suggest a fairly linear dependence of energy against probability of 
site occupation. 

For the first time in this work, it is necessary to explicitly average over the disorder, in this case the occupation 
numbers. This must be a quenched average and for this kind of site-disordered problem it is natural to use a formalism 
originally due to Morita [l5| and recently championed by Kuhn [1(1 . Essentially the method consists of taking an 
annealed average, but inserting constraints to fix the distribution of the occupation numbers at its quenched value. 
In practice this is achieved through an approximation in which only moments of the distribution are fixed. As a first 
step we would only consider the first two moments from the hierarchy corresponding to the constraints: 

M-x 

E n k = N (78) 

fe=0 

N 2 

2J n k n k , = d— (79) 

<kk'> 

The first constraint is straightforward and is incorporated simply by including the integration over the ao mode that 
was dropped in earlier sections since the set of n^'s chosen already obeyed the constraint. The second constraint 
involves a sum over nearest neighbours and prevents the occupied sites from clustering together thus minimizing path 
lengths. 

Z N (0)* E M E n ^-d— \ / ft d^e'S^E^ ( 8 °) 

n fc =0,l \<fefe'> / J k=0 p 

This approach relies on a better understanding of the unperturbed, N = M, lattice and takes us away from the aim 
of this paper, so we defer further study for the future. 
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FIG. 2: TSP on two sites. 



VI. CONCLUSIONS 

We have introduced a well denned discrete problem: the site disordered TSP with N cities on M sites, and have 
studied it using a saddle point approach. 

• In the coarse, N 3> M limit, we obtain a discrete version of the continuum theory we developed previously 
using a functional approach. The discrete approach shows that M saddle points contribute to the result rather 
than the single one identified in the functional work, but that (at least for constant density of cities) they 
are all equivalent. The effect of fluctuations can also be computed using the discrete approach. However, in 
the limit where there is a diagonal distance much shorter than neighbour distances and this correction might 
give information about the minimal path length in one dimensional problems, the saddle point approach suffers 
technical difficulties and gives incorrect results. The nature of the difficulties were exposed in the M = 2 example 
as being associated with additional saddle points becoming relevant and saddle points approaching branch cuts. 

• The saddle point approach to the ordinary N = M TSP has the correct scaling to be able to access the interesting 
low temperature regime even at the leading order. However, this fails to give correct results since the saddle is 
shallow and O(N) modes all contribute to the fluctuations about the saddle generating a correction that is of 
the same order as the leading term thus invalidating the approach. 

• We have discussed the dilute model and speculated on an approach based on an approximation originally due 
to Morita. 

• We expect that it would be possible to use the saddle point method to analyze models that are intermediate 
between the N ^> M regime and the N — M case with city densities that vary as N a for some power < a < 1 . 
The quadratic fluctuations in these models will be of order TV 1- " and generally allow an expansion in N~ a , but 
we anticipate the same problems with the low temperature limit. 

Of course a saddle point approach is not the only way to proceed and one interesting aspect of the formalism we have 
presented is that it describes the path weights in terms of particular terms in a generating function thus resembling 
approaches to other combinatorial optimization problems. This approach can certainly be used to enumerate for small 
N. 

Acknowledgment: D.L. would like to acknowledge a Discipline Hopping Award from the EPSRC during which 
much of this work was done. 



APPENDIX A: EXPLICIT CALCULATION FOR M = 2 

The general derivation of the coarse regime given in the text did not label the eigenvalues globally and concentrated 
on the real positive solution of the saddle point equation. Explicit computations for particular small values of M 
illuminate the role of other eigenvalues and saddle points. 

Consider the case M = 2 (figure ^ where a diagonal entry is allowed in the distance matrix as otherwise the path is 
completely determined. This model can be reinterpreted as a 1-dimensional Ising model on a ring with N sites subject 
to the constraint that states must have zero overall magnetization. The analysis of the Ising model provides additional 
intuition, but a constrained transfer matrix approach repeats essentially the same technique as in the text. Writing 
d = e~" n for the diagonal term and n — e _/3Jjv for the neighbour term, the parameters of the system are then d/n 
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FIG. 3: Complex plane for TSP on two sites. For p_ = and in the region where n < d, where it is cheaper to revisit a site 
than to move to the other site. The branch points (BP), cuts and two saddle points (SP) are shown. The angle t is given by 
sint = n/d. 



and p- — (po — pi)/2 using the terminology of section Hill The partition function is the sum of the contributions from 
each of the two eigenvalues Z = no\ni\(Z + + ZJ) where, 

Z± =l n dzo n dzi z ?° z ? A+ (Al) 
/ 2mz 2mz x 1 ± K ' 

and A± are the two eigenvalues of the transfer matrix, 

T = J_( d*_njz5zt\ (A2) 
z zx \ ny/z zi dz J 

In this form, one of the integration variables simply imposes the constraint n Q + ni = N that we will generally make 
implicit. The constraint is therefore redundant and if we proceeded in this way the two saddle point equations would 
be linearly dependent and only set zq and z\ proportional to each other. While this approach merely leads to a volume 
factor, it is clearer to remove the constraint from the start. We therefore insert the factor 2 / Ki8(zQZ\ — 1) (the 2m is 
the Jacobean factor) into the partition function to obtain a one dimensional contour integral: 

<]Z z n p-' 2 \1 (A3) 



2iriz 

Where A± are now the two eigenvalues of the transfer matrix, 



T =ft;j 



That is: 



A± = i ^(z + i)±^d 2 (z-i) 2 +4n 2 j (A5) 

In this form it is possible to keep track of the same eigenvalues, labeled by ± or their relative magnitude at z = 1, 
and not swap between them as z varies around the integration contour. Note that this has involved a choice of sign 
in moving between (|A2[) and (|A4[) . for example, consider z = zq = z\ = — 1 (obeying the constraint) in each version 
of the transfer matrix. The branch points and saddle points in the complex plane are shown in figure [3J For n < d it 
is possible to choose a contour on the unit circle, but for n > d, the contour must be deformed away from the branch 
cuts. 

For certain values of the parameters the integrals can be computed directly. For d — or n = 0, a contribution 
is only obtained for certain choices of no, n\. We find d = is only possible if tiq — ri\ — N/2 in which case 



Z = 2(N/2)l 2 n N . For n — then one of no or n\ must vanish. These restrictions follow from the integrals but match 
the expectations from the original geometric problem since the choice of connections must make it possible for paths 
to visit nodes the required number of times. For n — d, one eigenvalue vanishes and we obtain Z = Nld N irrespective 
of no, n\ since there are a total of ./V! paths each of equal weight. The limit of n <C d is interesting as the corrections 
to the leading term provide information about the real (not coarse) TSP. In this limit, the diagonal path is short and 
is always preferred over the neighbour path. By expanding we obtain Z = no\ni\Nn 2 d N ~ 2 indicating that there are 
N paths that only flip once between nodes and this matches the combinatorial approach. 

Now we take it that neither d nor n vanishes and use the saddle point approach. The stationarity equation is: 

- + ^ = ° (A6) 
z A± dz 

which can be solved to find two saddle points for each eigenvalue. These saddle points depend on the sign of p_, but 
without loss of generality we can always take the sign to be positive. 



_2 



1 + 2g 2 T Vg 2 +g 4 (A7) 



Z sp ^ • » 

2 2 

2 n P- 



■ 9 = W^T) (A8) 

The upper and lower signs correspond to the two different eigenvalues, while the two separate saddle points for each 
eigenvalue arise from the sign of z sp . 

In general the saddle points are located at different points depending on the eigenvalue in question. Now consider 
the symmetric case of p_ = where g 2 = 0, where the location of the saddle points at z sp = ±1, is the same for each 
eigenvalue. For Z + involving the eigenvalue A+, the value of the eigenvalue at the saddle points z sp — ±1 is (±d + n) 
with eigenvector (1, 1) at each saddle. For the other integral, Z_, involving eigenvalue A_, the value of the eigenvalue 
at the saddle points z sp = ±1 is (±d — n) with eigenvector (1, —1) at each saddle point. It is already apparent that 
for each integral, one of the saddle points dominates the other (since we assumed d ^ 0), for Z + the dominant saddle 
is at z sp = +1, whereas for Z_ it is z sp = — 1. 

These explicit computations match the expectations from the general development in the text. The dominant saddle 
points come from the real positive solution of equation (|17|) . and the subdominant saddle points arise from a distinct 
solution of the equation (|17[) . The eigenvector associated with the dominant saddle point depends on which integral, 
Z±, is considered. 

The asymptotic value of integrals for Z± can be obtained from the dominant saddle point in each case. By expanding 
about the dominant saddle we find that the contour should pass over it in the imaginary direction, and the contribution 
from fluctuations can be evaluated. Each of Z± yields the same contribution since N must be even for the symmetric 
case under consideration. Overall we obtain: 



z =w 2 > p vJ^ + <WH^j (A9) 

This expression reproduces the exact result at n = d, but the limits d — > or n — + do not reproduce the known 
results. In such limits the simple saddle point calculation fails for technical reasons. The subdominant saddles would 
be expected to contribute (in which case the contour over them should be in the real direction) and moreover the 
branch points can approach the saddle points or the cuts can force the contour to be deformed. A more careful 
approach via a uniform asymptotic expansion [12[ may give better results. 

In summary, the coarse TSP on two points allows us to carefully track each eigenvalue and find that multiple saddles 
can occur for each eigenvalue. One saddle point dominates each integral associated with a certain eigenvalue. In the 
full partition function each of these eigenvalues contribute equally The eigenvector corresponding to the dominant 
saddle depends on the eigenvalue being tracked. Finally, the limit of removing elements of the distance matrix by 
making them large can be delicate in the saddle point approach. In particular we are unable to recover the true TSP 
as a subleading contribution to the limit from the saddle point. 

Already for the M = 3 generalization, the coarse TSP on 3 sites, the algebra to allow an eigenvalue to be tracked 
becomes excessive even for d = and with the aid of a computer. Additional isolated solutions of the saddle point 
equations in the text for Sj can be found provided d ^ 0. 



[1] M. Mezard and G. Parisi, Europhys. Lett. 2, 913 (1986). 



1G 



[2] W. Krauth and M. Mezard, Europhys. Lett. 8, 213 (1989). 

[3] D.S. Dean, D. Lancaster and S.N. Majumdar, J. Stat. Mech. L01001 (2005). 

[4] D.S. Dean, D. Lancaster and S.N. Majumdar, Phys. Rev. E 72 026125 (2005). 

[5] D.S. Johnson and L.A. McGeoch, in G. Gutin and, A. P. Punnen (Eds), The Traveling Salesman Problem and its Variations 

(Combinatorial Optimization Series), Kluwer, Boston (2002). 
[6] M. Ghosh, S.S. Manna and B.K. Chakrabarti, J. Phys A 21, 1483 (1988); A. Chakraborti, and B.K. Chakrabarti, Eur. 

Phys.J. B 16, 677 (2000). 
[7] Philip J. Davis. Circulant Matrices. Chelsea Pub Co; 2nd edition (1994). 
[8] R.S. Garfinkel, Oper. Res. 25, 741 (1977). 

[9] B. D. McKay and N. C. Wormald, European J. Combin., 11 (1990) 565-580; B. D. McKay, Combinatorica, 10 (1990) 
367-377. 

[10] M. Mezard, G. Parisi, and A. Zee, Nucl. Phys. B559, 689 (1999); A. Zee and I. Affleck, J. Phys: Condens. Matter 12, 
8863 (2000). 

[11] R.M. Karp and J.M. Steele, in E.L. Lawler, J.K. Lenstra, A.H.G. Rinnooy Kan and D.B. Shmoys (Eds), The Traveling 

Salesman Problem, Wiley, (1985). 
[12] N.G. De Bruijn, Asymptotic Methods in Analysis North-Holland, Amsterdam (1958); F.W.J. Olver, Asymptotics and 

Special Functions Academic Press, New York (1974); R. Wong, Asymptotic Approximations of Integrals Academic Press, 

San Diego (1989). 

[13] J. Vannimenus and M. Mezard, J. Physique Lett. 45 L1145 (1984). 

[14] A. Barvinok, E. Kh. Gimadi and A.I. Serdyukov, in G Gutin and, A. P. Punnen (Eds), The Traveling Salesman Problem 

and its Variations (Combinatorial Optimization Series), Kluwer, Boston (2002). 
[15] T. Morita J. Math. Phys. 5, 1401, (1964). 
[16] R. Kiihn, Z. Phys. 100, 231 (1996). 



17 



